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ABSTRACT  ( Maximum  200  words) 

The  Phase-I  COIL  effort  has  focused  on  extending  AeroSoft's  GASP  and  SENSE  software  for 
COIL-laser  analysis  and  design.  Important  improvements  include:  modifications  to  the  species 
diffusion  model  to  account  for  effective  diffusion  coefficients  and  pressure-driven  diffusion  effects, 
modifications  to  account  for  laser  power  extraction,  the  addition  of  a  COIL  chemistry  model,  and  the 
development  of  a  strategy  for  implementing  a  coupled  flow-solver/power-extraction  sensitivity 
analysis  method.  Results  for  test  problems  involving  helium  and  iodine  injection  showed  significant 
effects  from  the  pressure-diffusion  terms.  Fully-coupled,  power-on  laser  simulations  exhibited  an 
instability  resulting  from  an  uncoupling  of  grid  points  caused  by  the  bi-linear  interpolation  method 
used  to  transfer  data  back  and  forth  between  the  optics  mesh  and  the  flow-solver  mesh.  This 
problem  occurred  when  there  was  a  large  disparity  in  the  number  of  grid  points  for  the  flow-solver 
and  optics  mesh  in  the  lasing  cavity.  Stable  power-on  results,  however,  could  be  obtained  for  well- 
matched  flow-solver/optics  grid  systems.  Future  plans  to  implement  a  conservative  interpolation 
scheme  should  eliminate  this  problem. 
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1  Introduction 


Chemical  lasers  represent  an  important  component  in  many  advanced  weapons  systems  of  interest  to 
the  Air  Force,  e.g.,  the  airborne  laser  program  (ABL).  The  chemical  oxygen-iodine  laser,  or  COIL, 
is  a  short  wavelength,  high-power  chemical  laser  that  operates  on  an  atomic  iodine  laser  transition. 
COIL  was  invented  in  1977  at  the  Air  Force  Weapons  Laboratory,  now  a  part  of  the  Air  Force 
Research  Laboratory  at  Kirtland  Air  Force  Base,  New  Mexico.  AFRL/DE  continues  to  develop 
COIL  technology  using  as  a  primary  test  apparatus  the  10  kilowatt-class  device  called  RADICL 
(Research  and  Development  Iodine  Chemical  Laser) 

Research  and  development  programs  for  COIL  typically  employ  numerical  analysis  techniques 
such  as  computational  fluid  dynamics  (CFD)  to  provide  insight  into  laser  performance  and  design. 
A  common  goal  shared  by  these  programs  is  the  need  to  more  thoroughly  understand  the  underlying 
physical  processes  of  chemical  lasers  so  that  performance  may  be  optimized.  Enhanced  CFD  analysis 
tools  can  play  an  essential  role  in  this  regard.  High-resolution  CFD  simulations  can  aid  in  interpreting 
'and  explaining  measured  data,  make  predictions  for  operating  conditions  beyond  the  test  database, 
'and  provide  accurate  estimates  of  laser  power  output.  Beyond  predicting  the  performance  of  a  given 
laser  configuration,  new  design  tools  are  required  to  aid  scientists  in  tuning,  scaling  and  shaping  a 
practical  system.  With  these  needs  in  mind,  AeroSoft’s  Phase-I  research  and  development  program 
has  focused  on  the  following  three  goals: 

1.  To  develop  accurate  and  efficient  CFD  analysis  tools  for  modeling  COIL  lasers. 
State-of-the-art  CFD  methods  will  be  extended  to  include  new  capabilities  and  physical  models 
important  in  COIL-type  laser  systems  such  as  the  RADICL. 

2.  To  develop  high-level  software  tools  tailored  for  the  design  and  optimization  of 
COIL  lasers. 

The  continuous  sensitivity  equation  method  (CSEM)  for  CFD-based  design  will  be  extended 
to  include  new  capabilities  and  physical  models  important  in  COIL-type  laser  systems. 

3.  To  develop  a  framework  for  coupling  discipline-specific  sensitivity  analyses. 

Work  specific  to  this  contract  will  lead  to  new  sensitivity-analysis  tools  for  the  coupled  fluid- 
dynamic/power-extraction  systems  in  COIL-laser  designs. 

The  Phase-I  COIL  effort  has  focused  on  extending  AeroSoft’s  QASV  and  SENSE  software  for 
COIL  lasers.  QASV  is  a  multi-zone,  upwind,  characteristic-based  code,  which  solves  the  integral 
form  of  the  Reynolds- Averaged  Navier-Stokes  equations.  QASV  models  the  fluid  continuum  with  a 
finite  number  of  control  volumes  and  incorporates  general  thermo-chemistry  modeling  with  an  ex¬ 
tensive  thermo-chemical  database.  The  latest  version,  QASV  v4,  has  more  efficient  time-integration 
schemes  including  implicit  treatment  of  zonal  boundaries.  QASV  v4  uses  a  distributed-memory 
parallel  architecture  in  order  to  exploit  a  wide  variety  of  platforms.  These  schemes  greatly  improve 
the  efficiency  of  the  code,  reducing  the  amount  of  time  needed  to  perform  a  computation.  AeroSoft’s 
SENSE  software  is  based  on  the  continuous  sensitivity  equation  method  and  provides  flow-field  sen¬ 
sitivities  to  a  wide  array  of  flow  and  geometric  design  variables.  SENSE  is  available  as  a  commercial 
package  and  has  been  used  in  a  variety  of  applications  of  aerospace  interest  [3,  4]. 
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The  modeling  effort  during  Phase-I  included  modifications  to  &  AST's  thermo-chemical  database, 
improvements  to  the  species  diffusion  model,  and  modifications  to  account  for  laser  power  extrac¬ 
tion.  The  thermo-chemical  database  manager  has  been  modified  to  allow  chemical  species  to  be 
present  at  multiple  energy  states,  e.g .,  I  ^Pi/2)  and  I  (2^3/2);  a  feature  characteristic  of  COIL 
chemistry  models.  A  20-reaction,  10-species  COIL  chemistry  model  has  also  been  added  to  the 
thermo-chemical  database.  The  species  diffusion  model  in  QASV  v4  has  been  extended  to  account 
for  effective  diffusion  coefficients  as  well  as  the  pressure-driven  diffusion  effects  that  are  important 
in  COIL  lasers.  The  enhanced  diffusion  model  has  also  been  added  to  SENSE.  The  diffusion  mod¬ 
els  in  both  QASV  v4  and  SENSE  have  been  tested  for  injection  problems  involving  helium  and 
iodine  at  conditions  similar  to  the  RADICL  injector.  The  results  show  significant  effects  from  the 
pressure-diffusion  terms.  Software  modifications  have  also  been  made  to  QASV  v4  to  account  for 
laser  power  extraction  in  the  optical  cavity.  This  task  involved  adding  lasing  source  terms  to  the 
iodine  species  continuity  equations  and  the  global  energy  equation,  implementing  a  geometric  optics 
resonator  model,  and  including  interpolation  methods  between  the  flow-solver  grid  and  the  optics 
grid.  Power-off  simulations  have  been  performed  for  the  RADICL  test  configuration  to  test  the  new 
chemistry  model.  Fully-coupled,  power-on  laser  simulations  exhibited  an  instability  resulting  from 
an  uncoupling  of  grid  points  caused  by  the  the  bi-linear  interpolation  method  used  to  transfer  data 
back  and  forth  between  the  optics  and  flow-solver  grids.  This  problem  occurred  when  there  was  a 
large  disparity  in  grid-point  density  between  the  flow-solver  and  optics  grids  in  the  lasing  cavity. 
Stable  power-on  results,  however,  could  be  obtained  for  well-matched  flow-solver/ optics  grid  sys¬ 
tems.  Future  plans  to  implement  a  conservative  interpolation  scheme  should  eliminate  this  problem. 
Phase-I  work  in  the  area  of  coupled  design  has  focused  on  defining  a  strategy  for  coupling  SENSE 
with  a  power-extraction  sensitivity  analysis  based  on  the  geometric  optics  resonator  model. 


2  Modifications  to  the  Thermo-Chemical  Database 

The  COIL  laser  operates  on  the  2(jPi/2)  -»  2  (P3/2)  transition  between  the  spin  orbit  components 
of  atomic  iodine.  This  transition  is  “pumped”  through  energy  transfer  from  the  excited  02(XA) 
state  of  molecular  oxygen.  Continuous  wave  operation  is  achieved  by  injecting  molecular  iodine  into 
a  primary  flow  containing  the  excited  states,  02  (L  A)  and  02(1£),  and  the  ground  state,  02 (3£), 
of  molecular  oxygen.  Once  the  molecular  iodine  has  been  injected  into  the  flow,  energy  pooling 
processes  involving  the  excited  states  of  oxygen  result  in  the  dissociation  of  molecular  iodine,  /2, 
into  atoms.  This  is  followed  by  rapid  energy  transfer  and  the  establishment  of  a  population  inversion 
in  atomic  iodine  according  to  the  “pumping”  reaction 

/  +  02(1A)->I* +02(3£),  (1) 

and  the  various  competing  mechanisms  that  work  to  reduce  I*.  Note  that  in  Eq  (1)  I  and  I*  denote 
the  2  (P3 /2 )  and  2(Pi/2)  states  of  atomic  iodine  respectively.  During  the  lasing  process,  it  is  the 
excited  I*  state  of  iodine  that  is  stimulated  to  emit  photons,  while  the  energy  in  this  system  is 
stored  in  the  02(*  A)  state.  Because  of  the  relatively  large  concentration  of  02(1A),  each  iodine 
atom  can  be  repumped  and  cycled  many  times  before  leaving  the  lasing  cavity. 

During  this  Phase-I  effort  the  COIL  chemistry  model  in  Table  1  has  been  added  to  Q AST's 
thermo-chemical  database.  The  /2  dissociation  mechanism  consists  of  the  two  separate  reaction 
paths  described  by  reactions  7  and  11.  Reactions  1,  2,  8,  9,  and  12  —  14  affect  the  J2  dissociation 
process  through  the  production  (or  destruction)  of  02(1S)  and  /2 .  Reaction  15  is  the  pumping 
reaction,  and  reactions  3  —  6,  10,  and  16  -  20  work  to  deactivate  the  system  by  reducing  either 
02 (x A)  or  /*.  The  forward  reaction  rates  for  each  reaction  are  constant  and  given  in  Table  1. 
The  backward  rates  are  computed  using  the  forward  rates  and  the  equilibrium  constant  which  is 
determine  through  the  LeRC  curve  fits  and  the  minimization  of  Gibb’s  free  energy. 

As  indicated  in  the  above  discussion,  the  COIL  chemistry  mechanism  is  used  (in  part)  to  model 
transitions  in  the  internal  energy  states  of  the  oxygen  and  iodine  species.  As  a  result,  the  various 
energy  states  of  a  given  molecule  or  atom  must  be  represented  as  a  separate  species  in  the  thermo¬ 
chemical  database.  To  accommodate  this  feature,  a  tagging  convention  has  been  added  to  the 
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No. 

Reaction 

Mechanism 

Kf(Kg  —  mole,mA ,  s) 

1 

CM1  A)  +  CM1  A)  ->  02(AE)  +  02(JE) 

I2  Dissociation  (+02(1S)) 

4.035  x  109 

2 

02(1S)  +  H2O  — >  O2OA)  +  H2O 

I2  Dissociation  (-02(1E)) 

9.635  x  109 

3 

CM1  A)  +  02(3E)  -t  02(3E)  +  o2(3e) 

Competing  (-02(1A)) 

2.409  x  103 

4 

02(1A)  +  H20  ->  02(3E)  +  H20 

Competing  (-02(1A)) 

9.635  x  102 

5 

02(1’A)  H-  CI2  — ^  02(^S)  -j-  CI2 

Competing  (-02(1  A)) 

4.818  x  10° 

6 

02(1  A)  +  He  -+  02(3S)  +  He 

Competing  (—  02(1A)) 

2.288  x  1010 

7 

I2  +  02{1  E)  ->  21  +  02(3E) 

I2  Dissociation  (+ 1) 

1.204  x  109 

8 

h  +  02(1S)  -»  J2  +  02(3E) 

I2  Dissociation  (“02(1E)) 

6.323  x  1010 

9 

h  +  02(1A)  - f  02(3£) 

I2  Dissociation  (+/2) 

4.215  x  106 

10 

j2 + /*  -*  j2* + / 

Competing  (—  /*) 

1.807  x  1011 

11 

/2*  +02CA)  -+  2  7  +  02(3E) 

I2  Dissociation  (+/) 

3.011  x  1010 

12 

72  +  02(3S)  — ►  72  +  02(3E) 

I2  Dissociation  (—  /|) 

1.807  x  1011 

13 

72*  +  7720  ->  72  +  H20 

I2  Dissociation  (-/2) 

2.409  x  109 

14 

72*  +  He  ->  I2  +  He 

/2  Dissociation  (-/2) 

6.022  x  105 

15 

/  +  02(1  A)  ->  J*  +  02(3E) 

Pumping  (+/*) 

6.624  x  107 

16 

J  +  OzCA)  ->  7  +  02(3E) 

Competing  (-02(1A)) 

9.635  x  106 

17 

7*  +  02(1A)  + 

Competing  (—  I*) 

1.626  x  104 

18 

7*  +  02(1A)  — >  7  +  CM1  A) 

Competing  (—  /*) 

4.697  x  1010 

19 

7*  +  /  ->  21 

Competing  (-/*) 

6.022  x  107 

20 

I*  +  H2O  — ^  1 H2O 

Competing  (—  /*) 

2.409  x  109 

Table  1:  10-species,  20-reaction  COIL  chemistry  model 


Species 

Database  Name 

02(3  E) 

0_2 

02{lA) 

(0_2)_[sd] 

02{l  E) 

(0_2)_ [ss] 

I2 

I_2 

4* 

1 — 1 

£ 

■P 

to 

1 _ 1 

1 

✓"~N 

CN 

1 

M 

I(2P3/  2) 

I 

I(2P1, 2) 

(I) _ [star] 

Table  2:  Thermo-chemical  database  species  naming  convention. 


thermo-chemical  database  manager  which  allows  a  chemical  species  to  be  present  at  multiple  energy 
states.  The  various  energy  levels  of  oxygen  and  iodine  have  been  entered  into  the  database  as  shown 
in  Table  2.  Figures  1  and  2  illustrate  the  new  naming  convention  as  shown  in  the  species  and  model 
windows  of  the  database  manager  for  the  COIL  mechanism. 


3  Effective  Diffusion  Model 

Pressure  diffusion  terms  have  been  shown  to  be  important  in  the  injector  region  of  the  COIL  where 
a  2:1  or  greater  pressure  ratio  must  exist  for  sonic  injection  to  be  achieved.  During  Phase-I  an 
effective  diffusion  model  with  pressure-gradient  terms  has  been  added  to  the  viscous  formulation 
in  QASV  v4.  The  diffusion  terms  appear  in  the  species  continuity  equations  as  well  as  the  global 
energy  equation.  The  continuity  equation  for  species  s  is  given  by 

it  III Padv+ £  (/°»v-“)dA=£(-^-“) dA+ Jff Usdv-  (2) 

where  ps  is  the  density  of  species  s,  V  is  the  velocity  vector,  u8  represents  the  chemistry  production 
term  for  species  s,  and  V,  is  the  diffusion  velocity  vector  for  species  s.  Prior  to  this  contract  the 
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Figure  1:  Species  window  illustrating  new  species  naming  convention  and  species  parsing 


Figure  2:  Model  window  illustrating  the  new  COIL  chemistry  model 
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(3) 


species  diffusion  velocity  was  given  by  Fick’s  law 

psVs  =  —pDsV(ps/p) , 

with  a  species  diffusion  coefficient  determined  using  a  constant  Schmidt  number  as 


Ds  = 


P 


pSc ' 


(4) 


The  effective  diffusion  model  represents  a  modeling  improvement  in  that  it  accounts  for  individual 
species  diffusion  coefficients  as  well  as  pressure-driven  diffusion  effects.  The  effective  diffusion  model 
for  species  s  is  given  by 


P«VS  —  —  Tt  Ms  D$m 


Vxs  + 


(*-7) 


Vp 
p  J 


N 


+  7 

H  j= 1 


jm 


,  (5) 


where  Ms  is  the  species  molecular  weight,  Dsm  is  the  species  effective  diffusion  coefficient,  is  the 
mixture  concentration,  and  Xs  is  the  species  mole  fraction.  Definitions  for  the  species  concentration, 
7S,  the  mixture  concentration,  and  the  species  mole  fraction  are  given  respectively  as 


Ps 

ls  =  77 
Ms 


N 

7 1  = 

$=  1 


7 1 


(6) 


The  species  effective  diffusion  coefficients  are  given  in  terms  of  a  mixture  summation  of  the  binary 
diffusion  coefficients,  T>ijy  as  follows 


DSm  — 


1  ~Xs 


Xj/Dsj 


and  satisfy  the  mass  conservation  constraint  that 


N 


'£psVs  =  o. 


(7) 


(8) 


S—l 


Closure  of  this  system  of  equations  requires  relations  for  the  binary  diffusion  coefficients,  Dy,  which 
have  been  determined  by  Hirschfelder  [5]  using  appropriate  collision  integrals.  The  resulting  equa¬ 
tions  are 


Vij  =  0.018828  T3/2 


(4-  +  -4-) 

\Mi  ^  Mj  ) 
pufMij 


(9) 


where 


<Ti  +  a  j 
=  — o —  > 


(10) 


a 


0.75571  0.68191 

+ 


2^*0.08742  1  y*0. 84910 


(11) 


and 


T 

T*  =  —  , 

Cij 


(12) 


(13) 
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Actual  implementation  of  the  effective  diffusion  model  has  been  performed  by  separating  the 
effects  of  mole-fraction  gradients  and  pressure  gradients.  This  feature  allows  the  option  to  run  with 
or  without  the  pressure  diffusion  terms.  Separating  these  effects  in  Eq  (5)  yields 

N 

PsV?  =  —  7t  Ms  Dsm  Vx«  +  7t  ^  E  Mi  **  >  (14) 

P  3= 1 

and 

PsV ?  =  — 7t (^-7)  ^ ^  (x* -  7)  y  -  (15> 

such  that 

P*VS  =  psV?  +  psVf  •  (16) 

Also  to  facilitate  ease  of  implementation,  the  following  variables  have  been  defined 

GX  =  MsD$mVXs,  (17) 


Gp  =  Ms  Dsm  Xs  - 


and 

Utilizing  these  definitions,  Eqs  (14)  and  (15)  become 

PsV?  =  -7 1 

and 

Ps^f  =  — 7t 


P* 


<55-7E<3? 

^  i=i 


TV 


G'-ef£Gt 

H  3=1 


Wp 

p 


(18) 


(19) 


(20) 


The  diffusion  term  resulting  from  mole-fraction  gradients,  Eq  (19),  behaves  in  a  similar  fashion 
to  Fick’s  law  and  results  in  molecular  diffusion  away  from  mole-fraction  gradients.  The  pressure 
diffusion  term,  Eq  (20),  results  in  molecular  diffusion  of  “lighter”  species  away  from  pressure  gradi¬ 
ents  and  “heavier”  species  toward  pressure  gradients.  To  see  this  we  consider  a  mixture  of  helium 
and  iodine  with  molecular  weights  M#e  =  4.003  and  M/2  =  253.82  respectively.  The  species 
densities,  pressure  and  temperature  are  taken  as  p#e  =  0.010395  Kg/m3,  pi2  =  0.007605  Kg/m?, 
p  -  9576  N/m2  and  T  =  438.5  K.  These  conditions  results  in  mass  fractions  of  PhJp  =  0.5775  and 
ph/p  =  0.4225  and  mole  fractions  of  XHe  =  0.988564  and  of  Xi2  =  0.011406.  Note  that  the  term 
(x»  -  P«/P)  in  Eq  (18)  will  be  positive  for  a  species  whose  mole  fraction  is  greater  than  its  mass 
fraction,  indicating  a  relatively  “light”  species.  This  same  term  will  be  negative  for  a  species  whose 
mole  fraction  is  less  than  its  mass  fraction,  indicating  a  relatively  “heavy”  species.  Evaluating  the 
terms  in  Eq  (20)  for  the  specified  conditions  yields 

PHeVpHe  =  -1.2456  x  10~4  ^  ,  (21) 


and 

phvvh  =  +1.2456  x  10~4  ^  .  (22) 

These  results  show  that  the  lighter  helium  molecules  diffuse  in  a  direction  opposite  the  pressure- 
gradient  vector,  while  the  heavier  iodine  molecules  diffuse  in  the  same  direction  as  the  pressure- 
gradient  vector.  The  net  effect  is  for  light  and  heavy  particles  to  tend  to  separate  in  regions  of  large 
pressure  gradients. 
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3-1  Helium-Iodine  Injection  into  an  Infinite  Chamber 

As  a  first  test  problem  for  the  new  diffusion  model,  we  considered  the  two-dimensional,  supersonic 
injection  of  helium  and  iodine  into  an  infinite  chamber.  Figure  3(a)  shows  the  81  x  161  computational 
mesh  used  for  this  case.  The  injector  face  extends  from  y  =  0.0  m  to  y  =  ±0.0005  m  with  flow 
entering  the  domain  from  left  to  right.  Solid  chamber  walls  extend  above  and  below  the  injector 
face.  Grid  points  are  clustered  in  the  horizontal  direction  near  the  wall  and  injector  face  and 
in  the  vertical  direction  near  the  injector  side  walls.  The  jet  Mach  number  was  M  =  1.3,  the 
temperature  was  T  =  288  K,  and  the  pressure  was  p  =  28728  N/m2.  The  species  densities  were 
pHe  =  0.047480  Kg/m3  and  pj2  =  0.034737  Kg/m?.  Initial  conditions  in  the  chamber  consisted 
of  pure  helium  at  Mach  number  of  M  =  0.01,  a  temperature  of  T  =  300  K,  and  a  pressure  of 
p  =  9576  N/m2.  The  initial  species  density  for  helium  was  pne  =  0.015368  Kg/m3.  The  jet 
conditions  were  imposed  using  a  split-flux  free-stream  boundary  condition  and  the  chamber  wall 
was  treated  using  a  no-slip  wall  boundary  condition.  All  other  boundaries  were  treated  using  first- 
order  extrapolation  from  the  interior. 

Figures  3(b),  3(c),  and  3(d)  show  the  resulting  Mach  number  contours  partially  superimposed 
with  streamtraces,  pressure  contours,  and  iodine  mass  fraction  contours.  Figures  3(b)  and  3(c) 
indicate  that  the  helium-iodine  mixture  undergoes  a  strong  expansion  as  it  leaves  the  injector, 
resulting  in  the  lower  pressures  and  higher  Mach  numbers  shown  by  the  contours  plots.  The  very 
strong  pressure  gradients  near  the  injector  side  walls,  i.e.,  y  =  ±0.0005  m,  result  in  large  pressure- 
diffusion  effects  which  work  to  separate  the  helium  from  the  iodine.  This  can  be  seen  in  Fig  3(d), 
where  the  iodine  concentrations  in  the  near-wall  region  are  essentially  zero.  Comparison  of  the 
streamlines  in  Fig  3(b)  and  the  iodine  mass  fraction  contours  in  Fig  3(d)  indicate  that  the  diffusion 
process  occurs  in  a  direction  transverse  to  the  mean  flow  velocity.  Figure  4  shows  the  normalized 
pressure  profile  and  helium  and  iodine  mass-fraction  profiles  for  the  vertical,  i  =  10  grid  line.  This 
grid  line  is  very  close  in  proximity  to  the  injector  face  and  chamber  walls.  The  normalized  pressure 
"profile  shows  a  large  pressure  gradient  at  y  =  ±.0005  m  near  the  injector  side  wall.  The  iodine  mass 
fraction  profile  shows  a  maxima  in  the  high  pressure-gradient  region,  while  the  helium  mass  fraction 
profile  shows  a  minima.  These  results  are  consistent  with  the  pressure-diffusion  effects  previously 
described.  As  helium  molecules  encounter  the  pressure  gradient  they  tend  to  diffuse  in  the  ±y 
direction  away  from  the  pressure  gradient,  resulting  in  a  minima.  Iodine  molecules  tend  to  diffuse 
in  the  —  y  direction  toward  the  pressure  gradient  resulting  in  the  observed  maxima. 

3.2  Transverse  Helium-Iodine  Injection  r 

As  a  second  test  problem  we  considered  the  two-dimensional,  transverse  injection  of  helium  and 
iodine  into  a  pure  helium  free  stream.  Figured  (a)  shows  the  113  x  81  computational  mesh  used 
for  this  case.  The  lower  surface  is  composed  of  a  solid  wall  with  the  injector  face  located  between 
x  =  0.0075  m  and  x  =  0.0077  m.  Grid  points  are  clustered  in  the  vertical  direction  near  the  wall 
and  injector  face  and  in  the  horizontal  direction  near  the  injector  side  walls.  The  jet  conditions 
are  identical  to  those  in  Section  3.1.  The  free-stream  conditions  axe  identical  to  the  initial  chamber 
conditions  in  Section  3.1,  except  the  Mach  number  is  M  =  0.8  with  the  flow  direction  being  left  to 
right.  As  before,  the  jet  conditions  were  imposed  using  a  split-flux  free-stream  boundary  condition 
and  the  lower  wall  was  treated  using  a  no-slip  wall  boundary  condition.  The  subsonic  inflow  boundary 
was  treated  by  fixing  the  all  the  flow  variables  except  pressure  to  the  free-stream  values.  The  pressure 
at  the  inflow  boundary  was  extrapolated  from  the  interior.  All  other  boundaries  were  treated  using 
first-order  extrapolation  from  the  interior. 

Figure  5(b)  shows  the  iodine  mass  fraction  contours  for  the  complete  domain.  Figure  6(a)  shows 
a  close-up  view  of  the  pressure  contours  and  streamtraces  in  the  injector  region.  Comparison  of  these 
figures  shows  three  distinct  flow  regions;  the  primary  flow  of  helium  which  expands  around  the  jet 
to  supersonic  speeds,  the  jet  flow  which  contains  regions  of  iodine  concentrations  that  are  slightly 
higher  than  jet  free-stream  values,  and  a  recirculation  region  just  to  the  right  of  the  injector  with 
lower  than  jet  free-stream  concentrations  of  iodine.  The  variation  in  iodine  concentrations  between 
the  jet  flow  region  and  the  recirculation  region  is  caused  by  the  pressure  diffusion  effects  that  tend 
to  separate  iodine  and  helium.  Figure  6(b)  shows  the  normalized  pressure  profile  and  iodine  mass 
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Figure  3:  Helium-Iodine  injection. 


Helium  and  Iodine  Mass  Fractions 

M=1 .3,  He-!2  Injection  Case  (1=10) 


Figure  4:  Helium  and  iodine  mass  fractions  (station  1=10) 


fraction  profiles  for  the  j  —  2  grid  line  in  the  injector  region.  The  normalized  pressure  profile 
shows  two  large  pressure-gradient  regions  in  the  vicinity  of  the  injector  side  walls.  The  iodine  mass 
fraction  profile  (for  the  effective  diffusion  model)  shows  peaks  at  these  locations,  indicating  that  the 
iodine  molecules  are  diffusing  in  the  direction  of  the  pressure-gradient.  This  causes  the  jet  flow  to 
liave  slightly  higher  concentrations  of  iodine  than  expected.  At  the  same  time  helium  molecules  are 
diffusing  away  from  the  pressure  gradient  causing  the  recirculation  region  to  be  composed  of  primarily 
helium.  Figure  6(b)  also  shows  the  iodine  mass  fraction  profile  for  a  simulation  using  Fick’s  law, 
which  has  no  pressure  terms.  In  this  case  both  the  jet  and  recirculation  region  have  helium  and 
iodine  concentrations  consistent  with  the  jet  free-stream  conditions.  Figures  7(a)  and  7(b)  show 
iodine  mass  fraction  contours  in  the  vicinity  of  the  injector  for  both  the  effective  diffusion  model 
and  Fick’s  law. 


4  Sensitivities  for  the  Effective  Diffusion  Model 

The  sensitivity  implementation  in  SZbfSZ  is  fundamentally  based  on  implicit-differentiation  applied 
to  the  boundary-value  problem  describing  the  fluid  mechanics.  That  is,  one  envisions  an  implicit 
equation 


ft(Q,fl=  0,  (23) 

where  Q  denotes  the  distributed  dependent  variables  ( e.g .  density,  momentum,  energy)  and  /? 
denotes  a  design  variable.  In  common  terminology  Q  is  known  as  the  state  and  p  as  the  design 
parameter.  For  fixed  /?,  Eq  (23)  is  solved  for  Q.  This  defines  a  map 

P  Q, 


which  associates  a  flow  solution  with  the  specified  design  parameter (s).  The  sensitivity  is  the  deriva¬ 
tive  of  this  map;  it  provides  a  linear  approximation  for  how  the  flow  solution  will  change  under  a 
small  change  in  the  design  parameter.  One  proceeds  by  formally  differentiating  Eq  (23)  to  produce 


dlZdQ  on 

dQ  dp  +  dp 


(24) 
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(a)  2-Dimensional,  113  x  81  grid. 


(b)  Iodine  mass  fraction  contours. 

Figure  5:  Helium-iodine  transverse  injection  case. 


The  S£Af$£  code  implements  a  numerical  approximation  to  Eq  (24)  for  situations  wherein  K  models 
Reynolds- Averaged  Navier-Stokes  flows  with  finite-rate  chemistry. 

The  integral  equations  for  the  fluid  dynamic  system  are  written  as 


Ftfffqdv  +  fA  (F(Q)  ■  fi)  dA  =  £  (F*(Q)  -VdA+ffJ  SdV>  (25) 


where  Q  =  Q(z,y,£,  t;  fi)  represents  the  vector  of  state  variables  resulting  from  conservation  of 
mass,  momentum  and  energy.  The  surface  integrals  represent  the  inviscid  and  viscous  fluxes  (F  and 
F v).  Changes  in  chemical  composition  are  governed  by  the  species  production  terms  in  the  source 
term  S.  Formal  differentiation  of  Eq  (25)  with  respect  to  the  generic  design  parameter,  /?,  results  in 
the  sensitivity  equations.  These  equations  are  linear  in  the  sensitivities.  Presented  in  integral  form, 
they  are 


FtfIICl'dV  +  fA(F'(Q)-h)dA  =  ^{F'v(Q)-n)dA  +  jjj  S '  dV, 


(26) 


where  Q'  is  the  flow  sensitivity,  i.e.,  dQ/dp.  The  inviscid  and  viscous  flux  sensitivities,  F'  and  F'r, 
and  the  source  term  sensitivity,  S',  are  formally  given  as 


dFdQ  dF^dQ  0S  5Q 

0Q  dp  ’  *  0Q  op  ’  dQ  dp  ' 


(27) 


S£MS£  solves  Eq  (26)  using  an  upwind  characteristic-based  formulation. 

As  a  first  step  toward  tailoring  S£MS£  for  COIL-laser  applications,  AeroSoft  has  extended  the 
viscous  sensitivity  formulation  to  include  the  effective  diffusion  model.  This  task  required  modifying 
the  species  diffusive  flux  sensitivity  terms  which  are  given  as 

K  =  -jp(p*Vs-n)  ,  (28) 
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(a)  Pressure  contours  and  streamtraces. 


Iodine  Mass  Fractions 

H e-l2  Transverse  Injection  Case  (J=2) 


(b)  Normalized  pressure  and  iodine  mass  fractions. 


Figure  6:  Helium-iodine  transverse  injection  case  (close-up  of  injector). 


(b)  Iodine  mass  fractions  for  Fick’s  law. 

Figure  7:  Helium-iodine  transverse  injection  case  (close-up  of  injector). 
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where  the  diffusion  velocity  is  determined  by  Eqs  (16)  (19)  and  (20).  The  sensitivity  formulation 
proceeds  by  formally  differentiating  Eqs  (19)  and  (20)  with  respect  to  the  design  variable,  yielding 


N 


g?“7  £g* 
H  j=  1 


(29) 


and 


N 


-it 


N 


G? 


7EGi 


where 

(psvs)' = (psv*y + (psv*)'  . 


(31) 


The  individual  sensitivity  terms  for  the  species  concentration,  the  mixture  concentration,  and  the 
species  mole  fraction  are  given  as 


N 


Ms 


5  =  1 


Xs  =  — 
It 


-  X.  -7 1 
Ms 


(32) 


,The  sensitivities  of  G*,  Gf,  and  Usro  are 

G*'  =  Ms  Z?;m  VXs  +  MsDsm^x's .  (33) 


and 


GPJ  —  Ms  D' 


+  Ms  Dsm 


D' 


Dsm 

1  —  X* 


N 

Xs  +  Dsm 

j=l, 


*1  T>'  . 
T>2  sj 


(34) 


(35) 


Finally,  closure  of  the  effective  diffusion  sensitivity  terms  is  given  by  sensitivities  of  the  binary 
diffusion  coefficient 


~  tJ  \2  T  Sly 


(36) 


where 


0.06606  0.57901 

rp  *1.08742  +  ji*1.84910 


n*t 


(37) 


and 

r'  =  —  .  (38) 

*U 


The  sensitivity  relations  given  by  Eqs  (29)  -  (38)  along  with  the  corresponding  flux-sensitivity 
Jacobian  terms  have  been  added  to  a  COIL  version  of  the  SENSE  software.  A  simplified  helium- 
iodine  chemistry  model  has  also  been  added  to  the  SENSE  package. 
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4.1  Sensitivities  of  Helium-Iodine  Injection 

The  sensitivity  analysis  has  been  applied  to  the  helium-iodine  injection  case  described  in  Section  3.1. 
For  this  case  the  design  variable  was  chosen  to  be  the  jet  iodine  density,  i.e.,  p  =  pi2Jet.  Fig¬ 
ures  8(a)  and  8(b)  show  the  iodine  mass  fractions  from  the  QASV  calculation  and  the  iodine  mass- 
fraction  sensitivities,  ( pi2/p )'  from  S£AfS£.  The  iodine  mass-fraction  sensitivities  at  the  injector 
face  are  imposed  by  the  jet  inflow  densities  and  density  sensitivities  through  the  following  relation 


The  jet  helium  and  mixture  densities  are  pi2  =  0.034737  Kg/m 3  and  p  =  .082217  Kg/m3.  The  jet 
helium  and  mixture  density  sensitivities  are  unity  for  this  case.  Under  these  conditions  Eq  (39)  gives 
a  value  of  ( ph/p )'  =  7.024;  the  free-stream  value  observed  in  Fig  8(b).  The  iodine  mass-fraction 
sensitivities  are  observed  to  have  values  larger  than  (pi2/p)'  =  7.024  near  the  injector  side-wall 
regions.  This  indicates  that  the  pressure-diffusion  effects  become  more  pronounced  as  the  iodine 
density  in  the  jet  increases.  To  see  that  this  is  plausible,  we  analyze  the  leading  pressure-gradient 
term  in  Eq  (15),  ( \h  ~  Ph/p)-  For  a  5%  increase  in  the  design  variable  the  absolute  value  of  this 
term  increases  from  |(x/2  -  Pi2 / p) I  =  0.387750  to  a  value  of  | (xi2  —  pi2/p) I  =  0.397994  (a  2.64% 
increase).  As  a  result  the  pressure  diffusion  effects  should  become  more  pronounced.  Figure  9(a) 
shows  iodine  mass  fraction  sensitivities  along  the  z  =  10  grid  line  for  both  S£MS£  and  a  central 
difference  approximation  using  QASV  and  a  5%  change  in  design  variable.  The  two  profiles  agree 
very  well  and  both  predict  an  increase  in  the  iodine  mass  fraction  near  injector  side  wall.  Figure  9(b) 
shows  a  close  up  of  iodine  mass  fraction  profiles  for  the  baseline  solution,  a  QASV  solution  for  a 
5%  increase  in  p,  and  a  Taylor’s  series  projection  using  the  baseline  flow  and  the  sensitivities  from 
S£MS£.  The  projected  solution  agrees  very  well  with  the  “exact”,  QASV  solution,  except  for  the 
^slight  lag  near  the  start  of  the  large  pressure-gradient  region. 

4.2  Sensitivities  of  Transverse  Helium-Iodine  Injection  /o 

l/o 

The  sensitivity  analysis  has  also  been  applied  to  the  transverse  helium-iodine  injection  case  described 
in  Sectioi/3.2.  For  this  case  the  design  variable  was  chosen  to  be  the  jet  velocity,  i.e.,  P  =  Viet. 
Figure fa)  shows  the  iodine  mass  fractions  from  the  QASV  calculation  and  Fig^qb)  shows  the  iodine 
mass-fraction  and  velocity-vector  sensitivities.  The  velocity- vector  sensitivities  indicate  an  increase 
in  the  jet  velocity  as  well  as  deeper  penetration  of  the  helium-iodine  mixture  into  the  primary  flow. 
The  iodine  mass-fraction  sensitivities  also  indicate  that  the  jet  flow  has  penetrated  further  into  the 
primary  flow.  The  region  of  positive  (and  negative)  sensitivities  just  above  (and  below)  the  jet 
indicate  that  the  jet  has  shifted  upward.  The  portion  of  the  recirculation  region  with  low  iodine 
mass  fractions  has  also  shifted  upward,  as  indicated  by  the  positive  sensitivities  in  this  area. 


5  Coupling  QASV  with  the  Power  Extraction  Model 

Current  COIL-CFD  analysis  codes  at  AFRL/DE  can  be  viewed  as  coupled  fluid  mechanics  and  laser 
power  extraction  codes.  The  MINT  implementation  [1]  of  the  COIL  simulation  couples  a  pde-based 
model  of  the  mechanics  for  compressible,  viscous,  chemically  reacting  flow  with  a  ray-tracing  model 
[2]  for  the  power  extraction.  The  QASV  implementation  described  in  this  section  follows  the  same 
construct. 

5.1  Lasing  Source  Terms 

Power  extraction  in  the  lasing  cavity  manifests  itself  in  the  form  of  source  terms  in  the  iodine  species 
continuity  equations,  and  in  the  global  energy  equation.  The  species  continuity  equations  are  affected 
because  components  at  different  energy  states  are  represented  as  separate  species.  In  a  COIL  laser, 
iodine  in  the  excited,  /(2P1/2),  energy  state  is  stimulated  to  emit  photons  (i.e.,  stimulated  emission). 
Upon  emission  of  the  photon,  the  iodine  atom  assumes  the  lower,  I  i^Pz/z) ,  energy  state.  As  a  result, 
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(a)  Iodine  mass  fractions. 


(b)  Iodine  mass  fraction  sensitivities. 


Figure  8:  Sensitivities  for  Helium-iodine  injection. 


power  extraction  affects  the  concentration  of  these  two  states.  Power  extraction  also  accounts  for  a 
net  energy  loss  (through  emission  of  photons)  to  the  flow.  The  modified  species  continuity  equation 
is  given  as 

JJJ  psdV  +  J^{psV  ■  n)  dA  =  J^(-psVs  ■ n)  dA  +  JJJ  usdV  ±  JJJ  ^J^-dV,  (40) 

where  the  last  term  on  the  right-hand  side  of  Eq  (40)  is  the  power-extraction  source  term.  The 
source  term  is  composed  of  the  gain,  a,  the  two-way  average  intensity,  I,  the  molecular  weight,  Ms, 
and  the  energy  of  a  photon,  hu  =  90956  J/g  -  mole.  The  field  variable,  a,  represents  the  gain  and 
is  defined  as 

(4l) 

where 
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Y(m) 


Iodine  Mass  Fraction  Sensitivities 


(P=P,a4' 


(pij/p)' 


(a)  Iodine  mass  fraction  sensitivities  (station  1=10). 


Iodine  Mass  Fractions  (Taylor  Projection) 


Pu/p 

(b)  Projected  iodine  mass  fractions  for  5%  increase  in  ft  (Close 
up). 

Figure  9:  Sensitivities  and  Taylor  projection  for  helium-iodine  injection. 
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(b)  Iodine  mass-fraction  and  velocity-vector  sensitivities. 

Figure  10;  Sensitivities  for  transverse  helium  iodine  injection. 
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Here  iV/.  and  Ni  represent  number  densities  for  the  excited  and  ground  states  of  atomic  iodine 
and  (j)(v )  is  the  Voight  line-shape  function.  The  constants  A  and  A  are  taken  as  A  —  5.0  sec-1  and 
A  =  1.31527  x  10~4  cm.  The  plus  sign  in  Eq  (40)  is  used  for  the  ground  state  of  iodine  and  the  minus 
sign  is  used  for  the  excited  state.  No  other  species  in  the  model  will  have  a  source  term  contribution 
from  the  power  extraction  model.  The  global  energy  equation  will  have  a  similar  source  term,  -a  I, 
to  account  for  the  energy  removed  from  the  flow. 

The  two-way  average  intensity,  I,  is  provided  as  output  from  the  optics  resonator  model  (in  this 
case  the  ray-trace  model  of  Ref  [2]),  and  is  a  function  of  the  gain,  the  mirror  radii  and  reflectivities, 
the  geometry  of  the  optical  cavity,  and  various  transmission  and  diffractive  losses. 

5.2  The  Coupled  Model 

In  the  coupled  model  we  have  two  systems,  a  fluid-dynamic  system  denoted  by  Ilf  and  a  laser  power 
system  denoted  by  TZP 


Kf(uf,up)  =  0  (45) 

Kp(uf,up)  =  0.  (46) 

Here,  the  variable  uj  describes  the  fluid-dynamic  state  and  the  variable  up  describes  laser-power 
state.  An  overview  of  the  coupled  fluid-dynamic/power-extraction  analysis  procedure  is  as  follows: 

1.  The  CFD  flow  solver  uses  a  pseudo-time  iteration  to  drive  the  flow  residual  (i.e.  71  <■ )  to  zero. 
The  state-variable  from  the  laser-power  extraction  model,  generically  up,  represents  the  electric 
field  intensity,  I.  The  intensity  appears  in  the  fluid-dynamic  system  in  the  continuity  source 
terms  for  the  iodine  species  and  in  the  global  energy  equation  as  described  in  Section  5.1. 

-  2.  For  a  fixed  intensity  field  I  the  flow  iteration  proceeds  to  reduce  the  flow-residual.  The  flow- 

state  variable,  generically  Uf,  includes  the  flow  velocities,  thermodynamic  variables  and  species 
densities.  From  these  known  quantities  we  compute  a  gain  field,  a,  using  Eq  (41).  Values  in 
the  gain-field  are  computed  pointwise  from  the  flow-field  values. 

3.  The  gain-field  is  passed  to  the  power-extraction  module  for  the  calculation  of  an  updated 
intensity  field.  There  are  several  steps  in  this  linking  process  that  arise  from  modeling  choices. 
In  particular,  the  flow  simulation  is  based  on  a  unit-cell  approximation  and  imposes  a  periodic 
variation  of  the  flow- variables  in  the  direction  of  the  optical-axis  ( y  in  [2]).  Hence  the  gain  field 
is  averaged  in  the  direction  of  the  optical  axis  across  the  active  region.  The  gain  field  on  the 
flow-solver  mesh  is  then  interpolated  to  a  two-dimensional  polar  mesh  that  is  consistent  with 
the  ray-trace  algorithm.  At  present  a  bi-linear  interpolation  method  has  been  incorporated. 
The  interpolated  gain  is  then  passed  to  the  ray-trace  algorithm  in  the  array  storg(j  ,i). 

4.  For  the  given  gain-field  and  given  optical  cavity  geometry  there  is  not  necessarily  an  associated 
steady  intensity  field.  Steady-state  lasing  occurs  for  the  ray-trace  algorithm  when  the  round- 
trip  gain-equals-loss  condition  is  satisfied.  Such  a  steady  field  is  achieved  only  at  gain  saturation 
where  gain-production  and  laser  extraction  are  in  equilibrium.  As  a  result,  the  flow  and 
power  models  must  be  solved  simultaneously.  In  order  to  make  sense  of  a  sequential  iteration 
procedure  one  must  define  a  unique  intensity  field  that  is  associated  with  the  given  gain  field. 
The  power  extraction  module  based  on  [2]  does  this  by  calculating  a  return  intensity  and  using 
an  under-relaxation  update.  This  process  is  described  in  more  detail  in  Section  5.3. 

5.  The  updated  intensity  field  from  the  power  extraction  module  is  communicated  back  to  the 
flow  model  through  an  array  of  intensities  stori  ( j  ,  i) .  These  intensities  are  then  interpolated 
from  the  polar  optics  mesh  to  the  flow-solver  mesh  where  they  are  used  to  compute  updated 
lasing  source  terms. 


18 


5.3  Power  Extraction 


The  main  work  of  the  power-extraction  module  is  done  within  the  routine  stable. F  which  was 
supplied  by  AFRL/DE.  A  significant  code  fragment  is  reproduced  here. 

1000  continue 
c 

c  calculation  of  intensity  on  subsequent  iterations,  i.e.,  iopt  =  0 
c 

sdifmax  =0.0 
do  1500  i  =  l,nthetal 
do  1200  j  =  1 , nrmax 

rod(j,i)  =  RABS(rmaxls(i)  *  (xdist  -  xa(j))  /  distm) 
sum  =  storg(j,i)  -  scripl 
alphaq(l)  =  storg(j,i) 
do  1100  iq  =  2,m 

alam  =  SQRT(1.0  +  (rod(j,i)  *  afl(iq))**2) 
alphaq(iq)  =  xratio(j,iq)  *  (storg(iset (j , iq)-l , i)  - 
1  storg(iset(j ,iq) ,i))  +  storg(iset (j , iq) ,i) 

sum  =  sum  +  alam  *  (alphaq(iq)  -  scripl) 

1100  continue 

s(j)  =  2.0  *  glength  *  sum 

exps  =  EXP(RMIN(expmin,s(j))) 

c  ensure  that  stori(j,i)  is  not  changed  on  restart 

if(ifirst  .eq.  1)  then 
iret(j)  =  inot (j ,i) 

.  else 

iret(j)  =  inot(j,i)  *  rlr2m  *  exps  *  prodj(j) 
end  if 

inot(j,i)  =  coninot  *  iret(j)  +  (1.0  -  coninot)  *  inot(j,i) 
expl  =  EXP (RMIN(expmin, glength  *  (alphaq(l)  -  scripl))) 
terml  =  (expl  -  1.0)  /  (glength  *  (alphaq(l)  -  scripl)) 
term2  =  1.0  +  prodj(j)  *  (rlr2m  *  exps  /  (rml  *  expl))  / 

1  (1.0  -  del(j) ) 

stori(j,i)  =  inot(j,i)  *  terml  *  term2 
c 

(several  lines  omitted) 

1200  continue 
c 

(several  lines  omitted) 

1500  continue 

The  i  and  j  loops  are  stepped  through  the  gain  region  using  polar  coordinates  to  span  the 
(x,y)  -  plane  in  a  coordinate  system  that  is  convenient  for  the  geometry  of  the  optical  cavity.  Each 
(i ,  j)  pair  labels  a  point  on  the  out-coupling  mirror  and  the  calculations  within  the  loop  implement 
Eq  (7)  in  [2,  page  5].  A  ray  of  initial  intensity  Io  leaving  a  point  on  the  flat  out-coupling  mirror 
in  an  orthogonal  direction  will  return  (orthogonally)  to  the  same  point  after  2 m  passages  through 
the  gain  medium.  After  each  passage  there  is  a  specular  reflection  from  the  opposite  mirror.  The 
return  intensity  is  given  by: 


Ir(J )  =  Io(J)(RiR2)m  exp  s  J]  (1  -  Sp)  (47) 

P=  1 
m 

s  =  2Lg^P-C)X,  (48) 

p=i 
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where  i?i  and  i?2  are  the  reflectivities  of  mirror  1  and  2,  and  ap  is  the  average  gain  along  the  p-th 
ray  segment.  This  calculation  is  completed  in  the  if  construct,  that  begins  with  the  line 

if(ifirst  .eq.  1)  then 

The  initial  intensity  field,  Jo,  for  the  next  invocation  of  this  routine  is  computed  in  the  line: 

inot(j,i)  =  coninot  *  iret(j)  +  (1.0  -  coninot)  *  inot(j,i) 

Note  that  this  code  fragment  is  embedded  in  an  i-loop  so  the  i  variable  does  not  appear  explicitly 
in  iret().  In  the  present  version  we  have  coninot  =  0.5,  so  that  the  scheme  under-relaxes  the 
intensity  field.  In  addition,  there  is  a  possibility  of  invoking  the  stable .  F  routine  several  times  with 
the  gain-field  frozen  at  the  current  QASV supplied  distribution.  When  the  coupled  system  converges 
the  initial  and  return  intensity  values  will  agree,  so  that,  in  effect,  the  relaxation  calculation  is 
identical  to  Io  =  I r}  or 

inot(j,i)  =  iret(j) 

Finally,  the  average  two-way  intensity  in  this  region  of  the  gain  field  is  computed  from  Eq  (29) 
in  [2,  page  11]. 


h  = 


I0(J)[exp[(a(J)-£)Lg]-l] 


1  + 


exp  [-  (a(J)  -  £)  Lg\ 


(49) 


[a(J)-C]Lg  [  Ri[l-S(J)] 

This  calculation  is  (almost)  implemented  in  the  final  assignment: 
stori(j,i)  =  inot(j,i)  *  terml  *  term2 

The  stori  array  is  communicated  back  to  the  QASV  flow  code  where  it  is  projected  onto  the 
'flow-solver  grid. 


5.4  Software  Modifications 

During  Phase  I,  AeroSoft  has  performed  modifications  to  QASV  v4  to  facilitate  the  power-extraction 
model.  The  new  additions  to  QASV  include: 

1.  Memory  allocation  for  the  gain  and  intensity  fields. 

2.  Gain  computation  routine. 

3.  Lasing  source  term  and  Jacobian  routines. 

4.  Parallel  logic  to  gather  the  gain  from  zones  (within  the  lasing  cavity)  on  distributed  processors. 

5.  Parallel  logic  to  scatter  the  intensity  to  zones  (within  the  lasing  cavity)  on  distributed  proces¬ 
sors. 

6.  Routines  to  define  the  optics  parameters  and  generate  the  polar  mesh. 

7.  Implementation  of  the  stable. F  routine. 

8.  Interpolation  routines  between  the  flow-solver  mesh  and  the  polar  optics  mesh. 

9.  Iteration  strategy  for  updating  intensity. 

10.  Flow  visualization  capability  for  gain. 
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5.5  Fully-Coupled  COIL  Laser  Simulations 

Initial  attempts  to  perform  a  fully-coupled,  power-on,  RADICL  simulation  yielded  solutions  with 
divergent  lasing  power.  The  resulting  flow  fields  exhibited  localized  regions  within  the  lasing  cavity 
where  the  intensity  continued  to  grow  unbounded.  It  was  hypothesized  that  this  instability  resulted 
from  an  “uncoupling”  of  grid  points  caused  by  the  bi-linear  interpolation  method  used  to  transfer 
the  gain  and  intensity  back  and  forth  between  the  flow-solver  mesh  and  the  polar  optics  mesh. 
Figure  11  shows  the  superimposed  flow-solver  and  polar  optics  meshes  and  will  be  used  to  illustrate 
the  potential  problems  with  the  bi-linear  interpolation  method.  Notice  for  the  particular  case  shown 
here,  the  flow-solver  mesh  is  very  coarse  relative  to  the  optics  mesh.  As  an  example,  we  consider  the 
step  of  transferring  the  intensity  from  the  optics  mesh  points  to  the  flow-solver  cell  center  (denoted 
by  the  solid  black  circle).  In  order  the  interpolate  the  intensity  from  the  polar  mesh  to  the  flow-solver 
cell  center,  the  bi-linear  scheme  finds  the  four  closest  surrounding  cells  in  the  optics  mesh  (he.,  the 
grey  squares).  These  four  points  are  then  used  in  a  bi-linear  interpolation  to  determine  the  cell- 
center  value  for  the  intensity.  However,  for  the  case  shown  in  Figure  11  there  are  two  mesh  points 
(the  red  triangles)  inside  the  flow-solver  cell  that  are  never  used  in  the  interpolation  procedure.  The 
intensity  at  these  “uncoupled”  points  has  no  influence  on  the  flow-field.  Points  similar  to  these,  in 
our  initial  RADICL  simulations,  coincided  with  locations  where  the  intensity  was  allowed  to  grow 
unchecked. 

In  order  to  analyze  this  potential  problem  with  the  coupled  model,  a  simplified  test  case  was 
devised.  In  developing  the  grids  for  this  case,  care  was  taken  to  distribute  the  grid  points  such  that 
a  nearly  one-to-one  mapping  between  flow-solver  and  optics  mesh  points  was  obtained.  By  keeping 
the  grid  point  distribution  nearly  the  same  for  the  two  meshes,  the  uncoupling  problem  should  be 
eliminated.  The  simplified  test  case  is  described  in  the  following  section. 


Figure  11:  Bi-linear  interpolation  scheme. 


5.5.1  Simple  Test  Case 

Power-on  simulations  have  been  performed  for  a  simplified  flow  geometry,  but  with  flow  conditions 
and  optical  parameters  consistent  with  the  RADICL  laser.  This  test  problem  reduced  the  debugging 
effort  and  served  to  demonstrate  the  coupled  processes  involved  in  the  new  power-extraction  model. 
Figure  12(a)  shows  the  grid  topology  for  the  151  x  51  flow-solver  mesh  and  the  101  x  51  polar 
optics  mesh.  The  number  of  grid  points  shown  in  this  figure  has  been  reduced  to  facilitate  viewing. 
Notice  that  the  x-direction  stretching  at  the  leading  and  trailing  edges  of  the  optics  mesh  is  nearly 
identical  to  the  flow  solver  mesh.  The  start  of  the  aperture  is  at  x  -  .0 6  m  and  the  aperture 
length  is  La  =  ,05  m.  The  aperture  height  is  Ha  -  .028  m.  The  distance  between  the  curved 
mirror  and  the  flat  out-coupling  mirror  is  Dm  =  0.8  m  and  the  radius  of  curvature  of  the  curved 
mirror  is  r2  =  10.0  m.  The  coefficients  of  reflectivity  for  the  mirrors  are  Ri  -  0.90  and  R2  =  0.99 
for  the  out-coupler  and  curved  mirrors  respectively.  The  gain  length  was  taken  as  Lg  =  0.25  m. 
Distributed  losses  in  the  gain  medium  and  diffractive  losses  at  the  aperture  were  neglected.  The 
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Variable 

High-Gain  Region 

Low-Gain  Region 

PHe 

0.00260088  ~Kg]rr? 

0.000860848  ~Kgjm? 

Po2 

0.00208539  Kg/m3 

0.000811512  Kg/m3 

Po2(1a) 

0.00161352  Kg/m3 

0.000831930  Kg/m3 

Po2(3e) 

2.29196  x  10“6  Kg/m3 

8.38014  x  10"8  Kg/m3 

Ph2o 

0.000247298  Kg/m3 

0.000110298  Kg/m3 

Pi 2 

0.000235656  Kg/m3 

7.09078  x  10“5  Kg/m3 

Pq 

8.97101  x  10"7  Kg/m3 

8.90463  x  10~9  Kg/m3 

Pi 

4.74918  x  10“5  Kg/m3 

I  1.79366  x  10~6  Kg/m3 

Pi * 

0.000319047  Kg/m3 

4.33212  x  10-6  Kg/m3 

Pci2 

0.00122394  Kg/m3 

0.000545616  Kg/m3 

u 

881.8  m/s 

990.9  m/s 

V 

0.0  m/s 

0.0  m/s 

w 

0.0  m/s 

0.0  m/s 

p 

1124.6  N/m2 

1200.0  N/m2 

Table  3:  Inflow  conditions  for  laser  test  case. 


inflow  boundary  was  split  into  two  portions;  with  the  lower  stream  having  a  high  concentration  of 
I*  (high-gain  region),  and  an  upper  stream  with  very  little  I*  and  02(1A)  (low-gain  region).  The 
low-gain  region  was  designed  to  account  for  the  region  of  low  gain  observed  near  the  upper  wall  in 
the  power-off  RADICL  simulations.  The  inflow  conditions  for  the  test  case  were  chosen  from  actual 
conditions  (upstream  of  the  lasing  cavity)  for  the  power-off  RADICL  simulation,  and  are  listed  in 
Table  3. 

For  this  test  case  the  uncoupling  problem  did  not  occur  and  the  power-on  solution  produced  a 
-converged  value  for  the  lasing  power.  The  solid  line  in  Figure  13  shows  the  convergence  history 
for  the  lasing  power  in  this  case.  Figures  12(b)  and  12(c)  show  the  gain  and  intensity  contours 
for  this  simulation.  Values  for  the  inflow  gain  in  the  high  and  low-gain  regions  are  approximately 
a  =  1.065  1/m  and  a  =  0.00751  1/m  respectively.  The  maximum  gain  occurs  in  the  high-gain 
region  just  before  the  start  of  the  aperture  and  has  a  value  of  amax  —  1.308  1/m.  Past  the  start 
of  the  aperture,  the  gain  is  dramatically  reduced  through  interaction  with  the  lasing  source  terms. 
The  maximum  intensity  occurs  at  the  aperture  leading  and  trailing  edges  with  a  value  of  Imax  = 
3.6  x  109  W/m2.  The  power  output  for  this  case  converged  after  approximately  1000  iterations  and 
was  P  =  13.1  KW. 

In  a  second  test,  the  same  case  was  run  but  with  a  coarsened  flow-solver  mesh.  The  coarse 
flow-solver  mesh  was  obtained  by  sequencing  five  levels  in  both  the  x  and  y  directions,  yielding  a 
31  x  11  grid.  The  101  x  51  optics  mesh  from  the  previous  test  above  was  used.  The  grid  mismatch 
in  this  case  should  increase  the  possibility  that  grid-point  uncoupling  similar  that  shown  in  Fig  11 
will  occur.  As  expected,  the  lasing  power  for  this  case  did  not  converge.  The  dashed  line  in  Fig  13 
shows  the  lasing  power  history  for  this  case. 

One  possible  fix  to  the  uncoupling  problem  would  be  to  use  a  conservative  interpolation  method. 
In  this  scheme,  a  pseudo-cell  would  be  constructed  around  each  optics  grid  point.  The  intensity 
at  the  flow-solver  cell  center  would  then  be  influenced  by  each  optics  cell  that  overlaps  with  the 
flow-solver  cell.  Using  this  method,  all  optics  mesh  points  would  be  included  in  the  interpolation 
stencil. 

5.5.2  Minimum- Volume  RADICL  Simulation 

Based  on  the  known  coupling  problem  with  the  bi-linear  interpolation  method,  new  RADICL  simu¬ 
lations  have  been  performed  for  better-matched  flow-solver /optics  grid  systems.  Figure  14(a)  shows 
the  grid  topology  for  the  201  x  26  x  4  flow-solver  mesh  and  the  61  x  15  polar  optics  mesh.  As  in 
the  previous  example,  the  aperture  length  is  La  =  .05  m.  The  aperture  height  is  Ha  =  .03  m,  the 
distance  between  the  curved  mirror  and  the  flat  out-coupling  mirror  is  Dm  =  0.8  m  and  the  radius 
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(a)  Coarsened  flow-solver  mesh  and  polar  optics  mesh. 


(b)  Gain  contours  (1/m). 


(c)  Intensity  contours  ( Watts fm 2). 


Figure  12:  Test  laser  configuration  and  results. 
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Figure  13:  Lasing  power  vs.  iteration  number. 


of  curvature  of  the  curved  mirror  is  r2  =  10.0  m.  The  coefficients  of  reflectivity  for  the  mirrors  are 
Ri  =  0.90  and  R2  —  0.99  for  the  out-coupler  and  curved  mirrors  respectively.  The  gain  length  was 
taken  as  Lg  =  0.25  m.  Distributed  losses  in  the  gain  medium  and  diffractive  losses  at  the  aperture 
were  neglected. 

Figure  14(b)  shows  the  gain  contours  for  the  power-off  case  while  Figure  14(c)  shows  the  gain 
contours  for  the  power-on  case.  For  the  power-on  case,  the  maximum  gain  occurs  in  the  high-gain 
region  just  before  the  start  of  the  aperture.  Past  the  start  of  the  aperture,  the  gain  is  dramatically 
"reduced  through  interaction  with  the  lasing  source  terms.  The  maximum  intensity  occurs  at  the 
aperture  leading  and  trailing  edges  with  a  value  of  Imax  =  9.2  x  108  W/m2.  The  power  output  for 
this  case  was  approximately  P  =  8.64  KW.  Figure  15  shows  the  intensity  contours  in  the  lasing 
cavity. 

6  Coupled  Sensitivity  Analyses 

During  Phase-I  the  AeroSoft/ICAM  team  has  focused  on  defining  a  strategy  for  coupling  the  flow- 
field  sensitivity  solver  SSNS8  with  a  sensitivity  analysis  method  for  COIL  power-extraction.  The 
primary  objective  is  to  obtain  flow-field  and  lasing  sensitivities  for  the  coupled  model  described  in 
Section  5.  Current  COIL-CFD  analysis  codes  at  AFRL/DE  can  be  viewed  as  coupled  fluid  mechanics 
and  laser  power  extraction  codes.  The  QASV  implementation  described  in  Section  5  follows  this 
same  construct,  as  does  our  initial  study  of  the  sensitivity  calculations.  In  future  efforts  we  will 
study  a  continuous  model  of  the  power  extraction  based  on  a  paraxial  wave  equation  [7,  6]. 

While  it  is  possible  to  develop  and  implement  a  sensitivity-solver  specifically  for  the  COIL  system 
we  propose  to  study  a  general  scheme  for  coupling  sensitivity  analyses  for  two  distinct  disciplines, 
here:  fluid  mechanics  and  laser  power  extraction.  The  virtue  of  this  approach  is  that  it  leads  to 
re-usable  software,  and  can  be  carried  over  to  other  coupled  systems  of  interest  to  the  Air  Force. 

In  an  implicitly-coupled  model  we  have 

TZf(uf,up,  P)  —  0  (50) 

np(uf,up,l3)  =  0,  (51) 

where  the  iif  variable  describes  the  fluid  state,  the  up  variable  describes  laser-power  state,  and  P  is 
a  scalar  design  variable. 
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(a)  201  x  26  x  4  flow-solver  mesh  with  embedded  61  x  15  optics  mesh. 


(b)  Gain  contours  for  power-off  case. 


(c)  Gain  contours  for  power-on  case. 

Figure  14:  Minimum- Volume  RADICL  simulation  (power  on  and  off). 


The  corresponding  coupled  sensitivity  equation  is  described  by  the  system 


d  Kf  d  Tlf 
dUf  dup 

dnP  dUp 

a  up  j 


L  duf 


(52) 


One  approach  to  the  coupled  sensitivity  problem  builds  on  the  existing  SSAfSS  software  by  employ¬ 
ing  a  Jacobi  or  Gauss-Seidel  iterative  procedure  to  the  coupled  system  (52). 

An  alternative  formulation  begins  with  an  explicit  representation  of  the  laser  field  replacing  (51) 
with 


up  =  n(uf,p)  (53) 

Note  that  the  existing  power-extraction  module  can  be  viewed  as  a  method  for  (approximately) 
evaluating  the  function  H.  If  this  explicit  representation  is  substituted  into  (50)  we  obtain 


=  0. 

This  result  can,  in  turn,  be  implicitly  differentiated  to  yield  the  coupled  flow-sensitivity  equation 


'<9  71, 

dKf 

dW 

duf 

'dKf 

on 

571,1 

duf 

d  Up 

duf 

d/3 

3  tip 

'  dp 

dp  . 
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Figure  15:  Minimum- Volume  RADICL  intensity  contours  ( W/m 2). 

6,1  Sensitivity  Calculations 

In  Section  6  above  we  provided  two  forms  for  the  sensitivity  calculation:  Eq  (52)  based  on  a  general 
implicitly-coupled  formulation;  and,  Eq  (54)  for  a  partially  explicit  formulation.  The  structure 
suggested  by  the  stable. F  routine  is  an  admixture  of  these. 

Conceptually  there  are  several  distinct  intensity  fields  represented  in  stable. F. 

1.  an  initial  intensity  field  represented  by  the  array  inot, 

2.  a  return  intensity  field  represented  by  the  array  iret,  and 

3.  a  two-way  intensity  field  represented  by  the  array  stori. 

The  purpose  of  the  code  is  the  evaluation  of  the  two-way  intensity  field,  which  represents  the  coupling 
to  the  flow  code.  The  initial  and  return  intensities  are,  in  a  sense,  internal  field  variables  used  in  the 
calculation.  At  convergence  of  the  coupled  flow-intensity  systems  these  two  internal  fields  have  the 
same  value.  Specifically,  the  algorithm  computes  a  return  intensity-field  and  ultimately  a  two-way 
intensity  field  based  on  the  current  gain-field  and  an  initial  intensity-field.  The  latter  is  computed 
by  a  relaxation  method  from  the  previous  invocation  of  the  routine.  Thus,  we  suggest  the  abstract 
description 


<"' =«(«/,  «?,£).  (55) 

At  convergence  u£+1  =  We  compute  the  sensitivity  at  a  given  solution  for  the  coupled  flow- 
intensity  fields,  so  that  we  use  the  semi-explicit  form  of  Eq  (51) 


Up, /?)  —  Up  T~L{uf,Up^  /?). 

With  this  structure  the  coupled  sensitivity  Eq  (52)  can  be  written  as: 


dKf 

duf 
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o  Uf 


dKf 

d  Up 


d  uf 
dp 
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(56) 


(57) 


where  lp  is  the  identity  operator  on  the  state-space  where  up  lives.  The  form  (57)  is  the  approach 
we  are  taking  in  the  initial  COIL  sensitivity  studies. 
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For  now  we  assume  that  the  design  parameter  of  interest  (generically,  /?)  does  not  explicitly 
appear  in  the  power-extraction  module.  As  a  practical  matter  this  rules  out  cases  where,  for  example, 
the  geometry  in  the  optical  cavity  depends  on  /?.  Mathematically,  we  assume  that  the  function  Ti 
in  Eq  (56)  does  not  depend  explicitly  on  f3.  Our  plan  is  to  use  an  iterative  approach  to  solve  the 
coupled  linear  system  (57).  For  example,  a  Gauss-Seidel  type  implementation  would  be  of  the  form 


fdUA  +i  dUf 

\duf  )  f  d/3 


Sp+1 


d/3 
dU 


dj^lon 

duv  p 


duf 


cn+l  ,  ®  cn 

bf  +  9 


(58) 

(59) 


where  we  have  used  5/  and  Sp  to  denote  the  flow  and  laser-intensity  state  sensitivities,  respectively. 
That  is, 


s>  5 

q  =  ®  UP 

~  dp' 


Equation  (58)  will  be  solved  using  a  modification  of  AeroSoft/s  SENSE  code.  The  basic  flow 
model  describes  compressible,  laminar,  viscous,  chemically  reacting  flow.  For  the  COIL  application 
an  effective  diffusion  model  has  been  added  to  account  for  pressure  gradient  effects  (see  Section  4). 
The  COIL  chemistry  model  shown  in  Table  1  will  be  directly  coded  into  a  COIL-specific  version  of 
SENSE  in  the  near  future. 

The  final  feature  required  is  a  modification  to  include  the  term  Sp.  This  accounts  for  the 
effect  of  variations  of  the  laser  intensity  on  the  source  terms  in  the  flow  equations.  The  term  is  cal¬ 
culated  pointwise  throughout  the  flow  field  using  the  most  recently  calculated  intensity-sensitivity, 
which  is  communicated  to  the  flow-sensitivity  code  in  the  same  way  the  intensity  field  is  communi¬ 
cated  to  the  flow-solver. 

The  intensity-sensitivity  Sp  is  updated  as  in  Eq  (59).  The  procedure  for  implementing  this 
calculation  is  realized  in  a  new  piece  of  code,  based  on  the  stable. F.  At  an  abstract  level  the 
dependence  of  the  map  Ti  on  the  flow-field  is  naturally  decomposed  into  two  stages:  a  map  from 
flow-field  to  gain-field;  followed  by  a  map  from  the  gain-field  to  the  intensity-field.  In  symbols 


Ti(uf,up)  =  M(N(uf),up ), 

where  N  maps  the  flow-field,  pointwise,  to  the  associated  gain-field,  and  M  maps  the  gain-field, 
intensity  field  pair  to  the  associated  intensity-field.  The  required  Frechet  derivative  is  calculated  via 
the  chain-rule 

dU  _  d  Md  N 
duf  dug  duf 

The  N  map,  operating  pointwise,  is  an  ordinary  function  of  the  flow-field  data  (see  Eq  41). 
Flow-field  sensitivity  information,  such  as  sensitivity  of  the  various  number  densities  to  the  design 
parameter  are  available  from  the  the  modified  SENSE- code.  From  Eq  (41)  and  the  chain-rule 
(again)  we  can  explicitly  calculate  the  gain-field  sensitivity.  The  intermediate  matrix  in  this  chain- 
rule  is  our  calculation  of  the  term  As  in  the  coupled  flow-power  calculation  various  grid  systems 
must  be  reconciled.  Eventually,  we  compute  an  array  of  gain  sensitivities,  say  storg_sense  ([],[]). 
This  is  the  representation  of  the  sensitivity  of  the  gain-field  to  changes  in  the  parameter. 

The  function  M  maps  the  gain-field,  initial  intensity  field  to  an  new  intensity-field  and  is  imple¬ 
mented  in  the  stable. F  code.  This  is  represented  symbolically  in  Eq  (56).  Our  initial  implementa¬ 
tion  of  the  Gauss-Seidel  iteration  for  the  sensitivity  calculation  will  follow  the  approach  used  in  the 
coupled  flow  solution.  Specifically,  the  up  variable  in  Eq  (55)  is  identified  with  the  initial  intensity 
field,  which  at  convergence,  agrees  with  the  return  intensity  field.  Equation  (55)  is  implemented  in 
stable.F  in  the  line 
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iret(j)  =  inot(j,i)  *  rlr2m  *  exps  *  prodj(j) 


The  sensitivity  of  the  initial/retum  intensity  field  will  be  represented  on  the  cavity  grid  by  a  new 
array  inot .sense  ( [] ,  []).  The  implementation  of  the  second  Gauss-Seidel  update,  i.e.,  Eq  (59),  is 
symbolically  given  by: 

inot_sense([]  ,  [])  =  inot_sense ([],[] )  *  rlr2m  *  exps  *  prodj(j) 

+  inot(j,i)  *  [  rlr2m  *  exps  *  prodj(j)  ]_sense 

where  [  rlr2m  *  exps  *  prodj(j)  ]  .sense  denotes  the  /3-derivative  of  the  term  in  brackets. 

0 

—  [rlr2m  *  exps  *  prod j  ( j )] 

Under  our  simplifying  hypothesis,  that  $  does  not  explicitly  effect  the  cavity,  the  term  simplifies  to 

[r lr2m  *  prodj ( j)]  ^  [exps] 


In  these  expressions,  exps  is  the  exponential  term 

exp(s), 


where  the  argument  s  is  given  by  Eq  (47),  viz 


m 

s  —  2Lg  —  C)\. 

p=  i 


"The  /^-parameter  enters  here  through  the  gain-field,  that  is  in  the  variable  ap  which  is  computed 
from  the  stored  gain-field  storg.  The  required  expression  is 


9  [exp(s)] 

dp 


m 

exp(s)  Ap 

P=  1 


d  ap 

dp  * 


Since  ap  is  an  average  of  the  gain  along  a  particular  ray  through  the  gain-field  it  is  a  linear  function 
of  the  data  in  the  array  storg  ([],[]).  Its  ^-derivative  is  this  same  linear  function  of  the  entries  in 
the  gain-sensitivity  array  storg_sense(  [] ,  [] ) . 

The  two-way  intensity-field  is  computed  from  the  initial  intensity-field  according  to  (49).  Hence 
the  /3-derivative  of  the  two-way  intensity-field  can  be  computed  by  differentiating  the  expression  in 
(49).  The  result  will  involve  the  sensitivity  of  the  initial  intensity-field,  represented  by  the  array 
inot .sense  ([],[]),  as  well  as  the  explicit  gain-field  sensitivity,  represented  by  storg_sense  ([],[]). 
The  resulting  two-way  intensity-field  sensitivity  will  be  coded  in  a  new  array  (say)  stori_sense(  [] ,  [] ) 
This  is  communicated  back  to  the  flow  sensitivity  code  SENSE  in  the  same  way  that  the  two-way 
intensity-field  (stori(  [] ,  [] )  is  communicated  to  the  flow  solver,  QASV. 

We  will  modify  this  code  so  that,  when  an  appropriate  flag  is  set,  we  also  calculate  the  derivative. 
For  example,  one  fragment  of  the  code  is  (new  code  marked  by  asterisks): 

sum  =  storg(j,i)  -  scripl 
alphaq(l)  =  storg(j,i) 

*  if  (flag^sense)  then 

*  sum^sense  =  storg_sense(j ,i) 

*  alphaq_sense(l)  =  storg_sense(j , i) 

*  endif 

do  1100  iq  =  2,m 

alam  =  SQRT(1.0  +  (rod(j,i)  *  afl(iq))**2) 
alphaq(iq)  =  xratio(j,iq)  *  (storg(iset (j ,iq)-l,i)  - 
1  storg(iset ( j , iq)  , i) )  +  storg(iset (j , iq) ,i) 
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sum  =  sum  +  alam  *  (alphaq(iq)  -  scripl) 

*  if  (flag_sense)  then 

*  alphaq_sense(iq)  =  xratio(j,iq)  *  (storg_sense(iset(j ,iq)-l,i)  - 

*  1  storg_sense(iset(j ,iq) ,i))  +  storg_sense(iset (j , iq) , i) 

*  sum.sense  =  sum_sense  +  alam  *  alphaq_sense(iq) 

*  endif 

1100  continue 

We  continue,  supplying  additional  code  as  necessary  to  compute  the  array  stori_sense(  [] ,  [] ) , 
the  representation  of  the  sensitivity  of  the  intensity-field  to  changes  in  the  design  parameter. 
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